R Notebook

source('../settings/settings.R')
source('commonFunctions.R')

set.seed(43)
drive1 <- read.csv('../data/processed/analysis/TT1_Drive_1_PP.csv')
drive2 <- read.csv('../data/processed/Analysis/TT1_Drive_2_PP.csv')
drive3 <- read.csv('../data/processed/Analysis/TT1_Drive_3_PP.csv')
drive4 <- read.csv('../data/processed/Analysis/TT1_Drive_4_PP.csv', stringsAsFactors = T)

dfSeg <- data.frame(rep(1, nrow(drive4)), rep(2, nrow(drive4)), rep(3, nrow(drive4)), rep(4, nrow(drive4)))
names(dfSeg) <- c("Seg1", "Seg2", "Seg3", "Seg4")

combinedDf_Seg1 <- cbind(drive4, 
                    drive1$MeanPP_Seg0, 
                    drive2$MeanPP_Seg1, drive3$MeanPP_Seg1, 
                    drive2$MeanPP_Seg0, drive3$MeanPP_Seg0,
                    drive2$StdPP, drive3$StdPP,
                    dfSeg$Seg1
                  )
combinedDf_Seg2 <- cbind(drive4, 
                    drive1$MeanPP_Seg0, 
                    drive2$MeanPP_Seg2, drive3$MeanPP_Seg2, 
                    drive2$MeanPP_Seg0, drive3$MeanPP_Seg0,
                    drive2$StdPP, drive3$StdPP,
                    dfSeg$Seg2
                  )
combinedDf_Seg3 <- cbind(drive4, 
                    drive1$MeanPP_Seg0, 
                    drive2$MeanPP_Seg3, drive3$MeanPP_Seg3, 
                    drive2$MeanPP_Seg0, drive3$MeanPP_Seg0,
                    drive2$StdPP, drive3$StdPP,
                    dfSeg$Seg3
                  )
combinedDf_Seg4 <- cbind(drive4, 
                    drive1$MeanPP_Seg0, 
                    drive2$MeanPP_Seg4, drive3$MeanPP_Seg4, 
                    drive2$MeanPP_Seg0, drive3$MeanPP_Seg0,
                    drive2$StdPP, drive3$StdPP,
                    dfSeg$Seg4
                  )

names(combinedDf_Seg1) <- c(names(drive4), 
                       "PP_Dev_1_Turning",
                       "PP_Dev_2_Straight", "PP_Dev_3_Straight", 
                       "PP_Dev_2_Turning", "PP_Dev_3_Turning", 
                       "Std_PP_2", "Std_PP_3", "Segment")
names(combinedDf_Seg2) <- c(names(drive4), 
                       "PP_Dev_1_Turning",
                       "PP_Dev_2_Straight", "PP_Dev_3_Straight", 
                       "PP_Dev_2_Turning", "PP_Dev_3_Turning", 
                       "Std_PP_2", "Std_PP_3", "Segment")
names(combinedDf_Seg3) <- c(names(drive4), 
                       "PP_Dev_1_Turning",
                       "PP_Dev_2_Straight", "PP_Dev_3_Straight", 
                       "PP_Dev_2_Turning", "PP_Dev_3_Turning", 
                       "Std_PP_2", "Std_PP_3", "Segment")
names(combinedDf_Seg4) <- c(names(drive4), 
                       "PP_Dev_1_Turning",
                       "PP_Dev_2_Straight", "PP_Dev_3_Straight", 
                       "PP_Dev_2_Turning", "PP_Dev_3_Turning", 
                       "Std_PP_2", "Std_PP_3", "Segment")

combinedDf_Seg1$Subject <- paste0(as.factor(combinedDf_Seg1$Subject), ".S1")
combinedDf_Seg2$Subject <- paste0(as.factor(combinedDf_Seg2$Subject), ".S2")
combinedDf_Seg3$Subject <- paste0(as.factor(combinedDf_Seg3$Subject), ".S3")
combinedDf_Seg4$Subject <- paste0(as.factor(combinedDf_Seg4$Subject), ".S4")

combinedDf <- rbind(combinedDf_Seg1, combinedDf_Seg2, combinedDf_Seg3, combinedDf_Seg4)

combinedDf$Subject <- paste0("#", str_pad(combinedDf$Subject, 2, pad="0"))
combinedDf$Segment <- as.factor(combinedDf$Segment)

combinedDf_NoStressor <- combinedDf[combinedDf$Activity == "NO",]
combinedDf_Cognitive <- combinedDf[combinedDf$Activity == "C",]
combinedDf_Motoric <- combinedDf[combinedDf$Activity == "M",]

# combinedDf_NoStressor$Subject <- as.factor(combinedDf_NoStressor$Subject)
# combinedDf_Cognitive$Subject <- as.factor(combinedDf_Cognitive$Subject)
# combinedDf_Motoric$Subject <- as.factor(combinedDf_Motoric$Subject)

COLOR_NORMAL <- list(color='rgb(120,120,120)')
COLOR_COGNITIVE <- list(color='rgb(158,202,225)')
COLOR_MOTORIC <- list(color='rgb(58,200,225)')
COLOR_FAILURE <- list(color='red')

THRESHOLD_MILD = 0.07
THRESHOLD_EXTREME = 0.2

MARKER_LINE_MILD = list(color="blue")
MARKER_LINE_EXTREME = list(color="red")

yAxis <- list(title = 'Perinasal Perspiration (Log)', range=c(-0.3, 0.5))

fig_NoStressor <- plot_ly(combinedDf_NoStressor, x = ~Subject, y = ~PP_Dev_2_Straight, type = 'bar', name = 'Cognitive - Mean PP (Straight)', marker=COLOR_COGNITIVE) %>%
  add_trace(y = ~PP_Dev_3_Straight, name = 'Motoric - Mean PP (Straight)', marker=COLOR_MOTORIC) %>% 
  add_trace(y = ~PP_Dev_2_Turning, name = 'Cognitive - Mean PP (Turning)', marker=COLOR_COGNITIVE) %>% 
  add_trace(y = ~PP_Dev_3_Turning, name = 'Motoric - Mean PP (Turning)', marker=COLOR_MOTORIC) %>% 
  add_trace(y = ~PP_Dev, name = 'Failure - PP Deviation', marker=COLOR_FAILURE) %>% 
  add_segments(x="#1.S1", xend="#7.S4", y = THRESHOLD_MILD, yend = THRESHOLD_MILD, name="Threshold: Mild Change of PP",
                           line=list(color="blue", dash = 'dot')) %>%
  # add_segments(x="#01", xend="#41", y = THRESHOLD_EXTREME, yend = THRESHOLD_EXTREME, name="Threshold: Extreme Change of PP",
  #                          line=list(color="darkred", dash = 'dot')) %>%
  layout(yaxis = yAxis, barmode = 'group', title="Failure Driving \n Group=No Stressor")

htmltools::tagList(fig_NoStressor)

A marker object has been specified, but markers is not in the mode
Adding markers to the mode...
A marker object has been specified, but markers is not in the mode
Adding markers to the mode...

yAxis <- list(title = 'Perinasal Perspiration (Log)', range=c(-0.3, 0.3))

fig_Cognitive <- plot_ly(combinedDf_Cognitive, x = ~Subject, y = ~PP_Dev_2_Straight, type = 'bar', name = 'Cognitive - Mean PP (Straight)', marker=COLOR_COGNITIVE) %>%
  add_trace(y = ~PP_Dev_3_Straight, name = 'Motoric - Mean PP (Straight)', marker=COLOR_MOTORIC) %>% 
  add_trace(y = ~PP_Dev_2_Turning, name = 'Cognitive - Mean PP (Turning)', marker=COLOR_COGNITIVE) %>% 
  add_trace(y = ~PP_Dev_3_Turning, name = 'Motoric - Mean PP (Turning)', marker=COLOR_MOTORIC) %>% 
  add_trace(y = ~PP_Dev, name = 'Failure - PP Deviation', marker=COLOR_FAILURE) %>% 
  add_segments(x="#12.S1", xend="#3.S4", y = THRESHOLD_MILD, yend = THRESHOLD_MILD, name="Threshold: Mild Change of PP",
                           line=list(color="blue", dash = 'dot')) %>%
  # add_segments(x="#02", xend="#22", y = THRESHOLD_EXTREME, yend = THRESHOLD_EXTREME, name="Threshold: Extreme Change of PP",
  #                          line=list(color="darkred", dash = 'dot')) %>%
  layout(yaxis = yAxis, barmode = 'group', title="Failure Driving \n Group=Cognitive")

htmltools::tagList(fig_Cognitive)

A marker object has been specified, but markers is not in the mode
Adding markers to the mode...
A marker object has been specified, but markers is not in the mode
Adding markers to the mode...

yAxis <- list(title = 'Perinasal Perspiration (Log)', range=c(-0.3, 0.5))

fig_Motoric <- plot_ly(combinedDf_Motoric, x = ~Subject, y = ~PP_Dev_2_Straight, type = 'bar', name = 'Cognitive - Mean PP (Straight)', marker=COLOR_COGNITIVE) %>%
  add_trace(y = ~PP_Dev_3_Straight, name = 'Motoric - Mean PP (Straight)', marker=COLOR_MOTORIC) %>% 
  add_trace(y = ~PP_Dev_2_Turning, name = 'Cognitive - Mean PP (Turning)', marker=COLOR_COGNITIVE) %>% 
  add_trace(y = ~PP_Dev_3_Turning, name = 'Motoric - Mean PP (Turning)', marker=COLOR_MOTORIC) %>% 
  add_trace(y = ~PP_Dev, name = 'Failure - PP Deviation', marker=COLOR_FAILURE) %>% 
  add_segments(x="#24.S1", xend="#9.S4", y = THRESHOLD_MILD, yend = THRESHOLD_MILD, name="Threshold: Mild Change of PP",
                           line=list(color="blue", dash = 'dot')) %>%
  # add_segments(x="#05", xend="#31", y = THRESHOLD_EXTREME, yend = THRESHOLD_EXTREME, name="Threshold: Extreme Change of PP",
  #                          line=list(color="darkred", dash = 'dot')) %>%
  layout(yaxis = yAxis, barmode = 'group', title="Failure Driving \n Group: Motoric")

htmltools::tagList(fig_Motoric)

A marker object has been specified, but markers is not in the mode
Adding markers to the mode...
A marker object has been specified, but markers is not in the mode
Adding markers to the mode...

library(nlme)
combinedDf$Subject = as.factor(combinedDf$Subject)
combinedDf$Activity = as.factor(combinedDf$Activity)
combinedDf$PP_Dev_Group = ifelse(combinedDf$PP_Dev > THRESHOLD_MILD, 1, 0)

model = lme(PP_Dev ~ 
              abs(PP_Dev_2_Straight)
              + abs(PP_Dev_3_Straight)
              + abs(PP_Dev_2_Turning) 
              + abs(PP_Dev_3_Turning)
              + factor(Activity), 
            random=~1|Subject,
            data=combinedDf,
            method="REML")

# anova(model)
summary(model)

Linear mixed-effects model fit by REML
 Data: combinedDf 

Random effects:
 Formula: ~1 | Subject
        (Intercept)   Residual
StdDev:  0.06061321 0.02272995

Fixed effects: PP_Dev ~ abs(PP_Dev_2_Straight) + abs(PP_Dev_3_Straight) + abs(PP_Dev_2_Turning) +      abs(PP_Dev_3_Turning) + factor(Activity) 

 Correlation: 
                       (Intr) a(PP_D_2_S a(PP_D_3_S a(PP_D_2_T a(PP_D_3_T fc(A)M
abs(PP_Dev_2_Straight)  0.212                                                   
abs(PP_Dev_3_Straight) -0.017 -0.131                                            
abs(PP_Dev_2_Turning)  -0.142 -0.560     -0.333                                 
abs(PP_Dev_3_Turning)  -0.532 -0.125     -0.499     -0.012                      
factor(Activity)M      -0.498 -0.117      0.040     -0.167      0.151           
factor(Activity)NO     -0.356 -0.101      0.300     -0.292     -0.100      0.532

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-0.55213368 -0.36118802  0.01449327  0.22338495  0.72942632 

Number of Observations: 84
Number of Groups: 84 

plot(model)

model = lme(PP_Dev ~ 
              abs(PP_Dev_2_Turning)
              + factor(Activity), 
            random=~1|Subject,
            data=combinedDf,
            method="REML")

# anova(model)
summary(model)

Linear mixed-effects model fit by REML
 Data: combinedDf 

Random effects:
 Formula: ~1 | Subject
        (Intercept)   Residual
StdDev:  0.06060021 0.02272508

Fixed effects: PP_Dev ~ abs(PP_Dev_2_Turning) + factor(Activity) 

 Correlation: 
                      (Intr) a(PP_D fc(A)M
abs(PP_Dev_2_Turning) -0.460              
factor(Activity)M     -0.486 -0.261       
factor(Activity)NO    -0.421 -0.358  0.544

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-0.50898879 -0.36206730  0.01408431  0.24267504  0.66132090 

Number of Observations: 84
Number of Groups: 84 

plot(model)

model = lme(PP_Dev ~ 
              PP_Dev_2_Straight + 
              PP_Dev_3_Straight + 
              PP_Dev_1_Turning + 
              PP_Dev_2_Turning + 
              PP_Dev_3_Turning + 
              Std_PP_2 + 
              Std_PP_3 +
              factor(Activity), 
            random=~1|Subject,
            data=combinedDf,
            method="REML")

# anova(model)
summary(model)

Linear mixed-effects model fit by REML
 Data: combinedDf 

Random effects:
 Formula: ~1 | Subject
        (Intercept)   Residual
StdDev:  0.06098518 0.02286952

Fixed effects: PP_Dev ~ PP_Dev_2_Straight + PP_Dev_3_Straight + PP_Dev_1_Turning +      PP_Dev_2_Turning + PP_Dev_3_Turning + Std_PP_2 + Std_PP_3 +      factor(Activity) 

 Correlation: 
                   (Intr) PP_D_2_S PP_D_3_S PP_D_1 PP_D_2_T PP_D_3_T S_PP_2 S_PP_3 fc(A)M
PP_Dev_2_Straight  -0.204                                                                
PP_Dev_3_Straight  -0.058 -0.186                                                         
PP_Dev_1_Turning   -0.874  0.318    0.192                                                
PP_Dev_2_Turning    0.609 -0.751   -0.024   -0.646                                       
PP_Dev_3_Turning   -0.175  0.190   -0.662   -0.099 -0.319                                
Std_PP_2            0.458 -0.306   -0.149   -0.507  0.669   -0.282                       
Std_PP_3           -0.880  0.361    0.132    0.810 -0.781    0.192   -0.742              
factor(Activity)M  -0.095 -0.093    0.321    0.125 -0.021   -0.214    0.016 -0.079       
factor(Activity)NO -0.445 -0.027    0.278    0.344 -0.287   -0.072   -0.360  0.411  0.440

Standardized Within-Group Residuals:
          Min            Q1           Med            Q3           Max 
-0.6289604009 -0.2872063327  0.0004924481  0.2889330700  0.6496413406 

Number of Observations: 84
Number of Groups: 84 

plot(model)

Machine Learning

combinedDf$PP_Dev <- NULL

combinedDf$Subject <- NULL
combinedDf$Activity_NO <- ifelse(combinedDf$Activity == "NO", 1, 0)
combinedDf$Activity_C <- ifelse(combinedDf$Activity == "C", 1, 0)
combinedDf$Activity_M <- ifelse(combinedDf$Activity == "M", 1, 0)
combinedDf$Activity <- NULL

combinedDf$InSegment1 <- ifelse(combinedDf$Segment == 1, 1, 0)
combinedDf$InSegment2 <- ifelse(combinedDf$Segment == 2, 1, 0)
combinedDf$InSegment3 <- ifelse(combinedDf$Segment == 3, 1, 0)
combinedDf$InSegment4 <- ifelse(combinedDf$Segment == 4, 1, 0)
combinedDf$Segment <- NULL

combinedDf$Class <- ifelse(combinedDf$PP_Dev_Group == 1, T, F)
combinedDf$PP_Dev_Group <- NULL

# library(mefa)
# combinedDf <- rep(combinedDf, 10) 

n_folds <- 10
params <- param <- list(objective   = "binary:logistic", 
               booster          = "gbtree",
               eval_metric      = "auc",
               eta              = 0.1,
               max_depth        = 8,
               alpha            = 1,
               lambda           = 0,
               gamma            = 0.3,
               min_child_weight = 0.3,
               subsample        = 1,
               colsample_bytree = 0.5)
           
# XGBoost Model     
ml_data <- as.matrix(combinedDf %>% select(-Class))

xgb_m <- xgb.cv(   params               = param,
                  data = ml_data ,
                  label =  combinedDf$Class,
                  nrounds             = 500,
                  verbose             = F,
                  prediction          = T,
                  maximize            = T,
                  nfold               = n_folds,
                  metrics             = c("auc", "error"),
                  early_stopping_rounds = 100,
                  scale_pos_weight      = 3.05)

# xgb_m$evaluation_log[xgb_m$best_iteration,"test_auc_mean"]
xgb_m$evaluation_log[xgb_m$best_iteration,]

NA

library(pROC)

# it = which.max(xgb_m$evaluation_log$test_auc_mean)
# best.iter = xgb_m$evaluation_log$iter[it]
# best.iter 

plot(pROC::roc(response = ifelse(combinedDf$Class==T, 1, 0),
               predictor = xgb_m$pred,
               levels=c(0, 1)),
     lwd=1.5) 

Setting direction: controls < cases